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Abstract — Optimal design of water distribution networks, 
which are governed by a series of Unear and nonlinear equations, 
has been extensively studied in the past decades. Due to their 
NP-hardness, methods to solve the optimization problem have 
changed from traditional mathematical programming to modern 
intelligent optimization techniques. In this study, with respect to 
the model formulation, we have demonstrated that the network 
system can be reduced to the dimensionality of the number 
of closed simple loops or required independent paths, and 
the reduced nonlinear system can be solved efficiently by the 
Newton-Raphson method. Regarding the optimization technique, 
a discrete state transition algorithm (STA) is introduced to solve 
several cases of water distribution networks. In discrete STA, 
there exist four basic intelligent operators, namely, swap, shift, 
symmetry and substitute as well as the "risk and restore in 
probability" strategy. Firstly, we focus on a parametric study of 
the restore probability pi and risk probability p2. To effectively 
deal with the head pressure constraints, we then investigate 
the effect of penalty coefficient and search enforcement on the 
performance of the algorithm. Based on the experience gained 
from the training of the Two-Loop network problem, the discrete 
STA has successfully achieved the best known solutions for the 
Hanoi and New York problems. A detailed comparison of our 
results with those gained by other algorithms is also presented. 

Index Terms — Discrete state transition algorithm, water distri- 
bution network, intelligent optimization, NP-hardness. 

I. Introduction 

PIPES, hydraulic devices (pumps, valves, etc.) and reser- 
voirs are connected in a water distribution network in a 
complex manner. The physical behavior of a looped network is 
governed by a set of linear and nonlinear equations, including 
continuity and energy equations, and head loss functions. The 
overall planning tasks to be performed in water distribution 
networks consists of three kinds of problems: layout, design 
and operation. Although these problems are not independent 
with each other, they can be formulated and solved separately 
from a technical point of view since each one can be con- 
sidered as a parameter when others are being solved. In this 
work, we focus on the optimal design problem. 

Optimal selection of pipe diameters to constitute a water dis- 
tribution network respecting certain pressure requirements has 
been shown to be an NP-hard problem |1 1, mainly due to two 
reasons: nonlinear equations and discrete-valued diameters. 
A terribly clumsy method for designing pipe network is by 
enumeration or complete trial and error [2J. Traditional meth- 
ods are to linearize and relax the problem firstly to facilitate 

Xiaojun Zhou is witli tlie School of Science, Information Technology and 
Engineering, University of Ballarat, Victoria 3353, Australia and School of 
Information Science and Engineering, Central South University, Changsha 
410083, China (tiezhongyu2010@gmail.com) 

Manuscript received ** , 2013; revised **, 2013. 



the use of linear programming and nonlinear programming, 
and then they have to round off the solution to the nearest 
discrete diameters PI-IT). Such algorithms can not guarantee 
global optima and sometimes cause infeasible solutions. In the 
last few decades, intelligent optimization techniques including: 
genetic algorithm |[8)-pO), simulated anneaUng fTTl, shuffled 
complex evolution J 12), ant colony optimization jl3|, p4) , 
harmony search p3| , particle swarm optimization |[T6| , dif- 
ferential evolutioiipT) and some of their hybrids ]18| , p9) , 
have found wide applications in this field. The advantages of 
using these stochastic algorithms are: (1) simple representation 
of a discrete-valued solution; (2) independent to the problem 
structure to some extent; (3) easy computation due to the only 
use of the information about the objective function; (4) high 
probability to gain the global optimum or approximate global 
optimum in a reasonable amount of time. 

We introduce the recently developed intelligent optimization 
algorithm, state transition algorithm (STA) p0)-p2), which 
shows fantastic performance in continuous function optimiza- 
tion. In |J23), a discrete STA was proposed to solve the trav- 
eling salesman problem, and the results demonstrated that it 
consumed much less time and had better search ability than the 
well-known simulated annealing and ant colony optimization. 
The goal of this paper is to apply the discrete STA to the 
optimal design problem of the water distribution networks. 

This paper is organized as follows. In Section 11, the opti- 
mization model of water distribution networks is established, 
including the objective function, decision variables and some 
constraints. In Section III, the basic key elements in discrete 
STA are introduced. It focuses on the intelligent operators of 
discrete STA and a parametric study of the "restore proba- 
bility" and "risk probability" is the emphasis. How to deal 
with the constraints and the implementation of the discrete 
STA for the optimal design problem are illustrated in Section 
IV. In Section V, several case studies are given. The Two- 
Loop network is mainly studied to investigate the effect of 
penalty coefficient and search enforcement on the performance 
of the discrete STA. The gained experience is applied to 
other cases and the results achieved by the proposed discrete 
STA with other optimization algorithms are presented as well. 
Conclusion is derived in Section VI. 

II. Optimization model formulation of water 

DISTRIBUTION NETWORKS 

For a given layout of pipes and a set of specified demand 
patterns at the nodes, the optimal design of a water distribution 
network is to find the combination of commercial pipe sizes 
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which gives the minimum cost, subject to the following 
constraints: 

• continuity of flow; 

• head loss; 

• conservation of energy; 

• minimum pressure head. 



A. The objective function 

Considering that the pipe layout, connectivity and imposed 
minimum head constraints are known, in the optimal design 
problem of the water distribution network, the pipe diameters 
are the only decision variables. As a result, the objective 
function is assured to be a cost function of pipe diameters 



NP 



(1) 



where, f7 is a set of commercial pipe sizes, NP is the number 
of pipes, and Lj is the length of pipe j, which is known in 
this study. c{Dj) indicates that for every commercial pipe size, 
there is a corresponding cost per unit associated with it. 



B. Continuity equation 

Conservation of mass at nodes or junctions in a water 
distribution network yields a set of linear algebraic equations 
in terms of flows. At each node, flow continuity should be 
satisfied. 



^Q^n+Yl Qout +DM^O, 



(2) 



where, DM is the demand at the node, Qin and Qout are the 
flow entering and leaving the node, respectively. 



C. Head loss equation 

The head loss in a pipe in the water distribution network can 
be computed from a number of empirically obtained equations. 
The two commonly used equations are the Darcy-Weisbach 
head loss equation and the Hazen-Williams head loss equation. 
The general expression for the head loss in a pipe j located 
between nodes ; and k is given by 



Hi — Hk — Vj Qj I Qj 



Li 



C"D 



Qj I Qj 



(3) 



where. Hi and Hk are nodal pressure head at the end of the 
pipe at node / and k respectively; rj is called resistance factor 
for the pipe j; Qj is the flow in pipe j; cj is a numerical 
conversion constant depending on the units used; Lj is the 
length of pipe j; C is the roughness coefficient; a and /3 are 
coefficients. 

For International System of Units (SI), uj = 10.5088, a = 
1/0.54 = 1.85 and P = 2.63/0.54 = 4.87 are employed in 
this study using the Hazen-Williams formula. 



D. Energy equation 

Energy conservation equations around closed simple loops 
or between fixed head nodes along required independent paths 
in a network are nonlinear. Upon traversing a closed simple 
loop or a required independent path, the sum of pipe head 
losses around the loop or the path must be zero, which can be 
expressed as 

^' "i\Qor''-y,EL,^Q, (4) 



E 






where, Lg is the indices of pipes in a closed simple loop or a 
required independent path; ELj is the hydraulic grade line at 
the reservoir y. 

E. Minimum pressure head 

The minimum pressure head constraints at each node are 
given as follows 

ff^>if^mi„,V^=l,••• ,7VJ, (5) 

where. Hi ^nin is known, and NJ is the number of nodes. 

III. A BRIEF REVIEW OF THE DISCRETE STATE 
TRANSITION ALGORITHM 

Let's consider the following unconstrained integer optimiza- 
tion problem 



niin/(a;) 



(6) 



where, x = {xi, ■ ■ ■ ,x„), Xi E I C Z'"\ i = 1,- ■ ■ ,n, and 
f(x) is a real- valued function. 

A. The framework of the discrete state transition algorithm 

If a solution to a specific optimization problem is described 
as a state, then the transformation to update the solution 
becomes a state transition. Without loss of generality, the 
unified form of discrete state transition algorithm can be 
described as 

Xk+i = Akixk) Bk{uk) 
Vk+i = f{xk+i) 



(7) 



where, Xk e Z" stands for a current state, corresponding to a 
solution of a discrete optimization problem; Uk is a function 
of Xk and historical states; Ak{-), Bk{) are transformation 
operators, which are usually state transition matrixes; is a 
operation, which is admissible to operate on two states; / is 
the cost function or evaluation function. 

As a intelligent optimization algorithm, the discrete state 
transition algorithm have the following five key elements: 

(1) Representation of a solution. In discrete STA, we choose 
special representations, that is, the permutation of the set 
{1,2,--- ,ri}, which can be easily manipulated by some 
intelligent operators. The reason that we call the operators 
"intelligent" is due to their geometrical property (swap, shift, 
symmetry and substitute), and a intelligent operator has the 
same geometrical function for different representations. A big 
advantage of such representations and operators is that, after 
each state transformation, the newly created state is always 
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feasible, avoiding the trouble into rounding off a continuous 
solution in other cases. 

(2) Sampling in a candidate set. When a transformation 
operator is exerted on a current state, the next state is not 
deterministic, that is to say, there are possibly different choices 
for the next state. It is not difficult to imagine that all possible 
choices will constitute a candidate set, or a "neighborhood". 
Then we execute several times of transformation (called 
search enforcement SE) on current state, to sampling in the 
"neighborhood". Sampling is a very important factor in state 
transition algorithm, which can reduce the search space and 
avoid enumeration. 

(3) Local exploitation and global exploration. In optimiza- 
tion algorithms, it is quite significant to design good local 
and global operators. The local exploitation can guarantee 
high precision of a solution and convergent performance of 
a algorithm, and the global exploration can avoid getting 
trapped into local minima or prevent premature convergence. 
In discrete optimization, it is extremely difficult to define 
a "good" local optimal solution due to its dependence on 
a problem's structure, which leads to the same difficulty in 
the definition of local exploitation and global exploration. 
Anyway, in the discrete state transition algorithm, we define 
the slow change to current solution by a transformation as 
local exploitation, while the big change to current solution by 
a transformation as global exploration. 

(4) Self learning and regular communication. State transition 
algorithm behaves in two styles, one is individual-based, 
the other is population-based, which is certainly a extended 
version. The individual-based state transition algorithm fo- 
cuses on self learning, in other words, with emphasis on the 
operators' designing and dynamic adjustment (details given in 
the following). Undoubtedly, communication among different 
states is a promising strategy for state transition algorithm, as 
indicated in f22|. Through communication, states can share 
information and cooperate with each other. However, how to 
communicate and when to communicate are key issues. In 
continuous state transition algorithm, intermittent exchange 
strategy was proposed, which means that states communicate 
with each other at a certain frequency in a regular way. 

(5) Dynamic adjustment. It is a potentially useful strategy 
for state transition algorithm. In the iteration process of an 
intelligent algorithm, the fitness value can decrease sharply in 
the early stage, but it stagnates in the late stage, due to the 
static environment. As a result, some perturbation should be 
added to activate the environment. In fact, dynamic adjustment 
can be understood and implemented in various ways. For 
example, the alternative use of different local and global 
operators is dynamic adjustment to some extent. Then, we 
can change the search enforcement, vary the cost function, 
reduce the dimension, etc. Of course, "risk a bad solution in 
probability" is another dynamic adjustment, which is widely 
used in simulated annealing (SA). In SA, the Metropolis 
criterion |[24) is used to accept a bad solution: 



where, AE = f{xk+i) — f{xk), fcs is the Boltzmann 
probability factor, T is the temperature to regulate the process 
of annealing. In the early stage, temperature is high, and it 
has big probability to accept a bad solution, while in the late 
stage, temperature is low, and it has very small probability to 
accept a bad solution, which is the key point to guarantee the 
convergence. We can see that the Metropolis criterion has the 
ability to escape from local optimality, but on the other hand, 
it will miss some "good solutions" as well. 

In this study, we focus on the individual-based STA, and 
the main process of discrete STA is shown in the pseudocode 
as follows 
repeat 

[Best,fBest] <— swap(fcn,Best,fBest,SE,n,TOa) 
[Best,fBest] ^ shift(fcn,Best,fBest,SE,n,rnb) 
[Best,fBest] <— symmetry(fcn,Best,fBest,SE,n,mc) 
[Best,fBest] 4— substitute(fcn,Best,fBest,SE,set,n,mrf) 

> greedy criterion 



if fBest < fBest* then 

Best* <— Best 

fBest* ^ fBest 
end if 
if rand < pi then 

Best -s— Best* 

fBest -s- fBest* 
end if 
untU the maximum number of iterations is met 



> restore in probability 



As for detailed explanations, swap function in above pseu- 
docode is given as follows for example 



State -i— op_swap(Best,SE,n,r7ia) 
[newBest,fnewBest] -s— fitness(funfcn,State) 



if fnewBest < fBest then 
Best 4— newBest 
fBest ^ fnewBest 
else 

if rand < p2 then 
Best 4— newBest 
fBest -i— fnewBest 
end if 
end if 



[> greedy criterion 



> risk in probability 



probability p = exp( 



— ) 



(8) 



1 

2: 

3: 

4: 

5: 

6: 

7: 

8: 

9: 
10: 
11 

From the pseudocodes, we can find that in discrete STA, in 
the whole, "greedy criterion" is adopted to keep the incumbent 
"Best*", in the partial, a bad solution "Best" is accepted in 
each inner state transformation at a probability p2, and in 
the same while, the "Best*" is restored in the outer iterative 
process at another probability pi. The "risk a bad solution in 
probability" strategy aims to escape from local optimal, while 
the "greedy criterion" and "restore the incumbent best solution 
in probability" are to guarantee a good convergence. 



B. The representation, local and global operators 

In discrete STA, we use the index of the a commercial 
size as a representation for a solution to the optimal design 
problem. For example, if there are 8 pipes and for each pipe 
there are 3 choices, then the details of four special geometric 
operators are defined as follows 
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(1) Swap transformation 

Xk+l 



Al'"''^ima)xk, 



(9) 



where, A^"'' e Z"^" 



is called swap permutation matrix, ma 
is a constant integer called swap factor to control the maximum 
number of positions to be exchanged, while the positions 
are random. If m^ = 2, we call the swap operator local 
exploitation, and if nia > 3, the swap operator is regarded 
as global exploration. Fig. [T] gives the function of the swap 
transformation graphically when nia — 2. 



1 



12 3 3 13 



swap{2,6) 



I 3 



H 



1 



3 1 



2 I 



where, A|"^ g Z"^" is called substitute permutation matrix, 
md is a constant integer called substitute factor to control 
the maximum number of positions to be substituted. By the 
way, the positions are randomly created. If m^ = 1, we call 
the substitute operator local exploitation, and if nid > 2, 
the substitute operator is regarded as global exploration. Fig. 
|4] gives the function of the substitute transformation vividly 
when rriii = 1. 



l'l^|3|3|l|3|?1T 



substitute{2) 



i 



I I 1 I 3 I 3 I 1 I 3 I 2 I r 



Fig. 4. illustration of the substitute transformation 



Fig. 1. illustration of the swap transformation 



(2) Shift transformation 



* shift/ \ 

Xk+l = ^fc (mbjxk, 



(10) 



where, A 



shift 



e Z'' 



is called shift permutation matrix, rnf, 



is a constant integer called shift factor to control the maximum 
length of consecutive positions to be shifted. By the way, the 
selected position to be shifted after and positions to be shifted 
are chosen randomly. Similarly, shift transformation is called 
local exploitation and global exploration when nib = 1 and 
mb > 2 respectively. To make it more clearly, if nib = 1, we 
set position 2 to be shifted after position 6, as described in 
Fig. in 



i ,,,▼,,, '*'^?yf(2,6) 

I I 2 I 3 I 3 I I I 3 I 2 I l] ^ 



Fig. 2. illustration of the shift transformation 



(3) Symmetry transformation 



I I 3 I 3 I I I 3 I 2 I 2 I I 



Xk+l 

where, Af " e Z"><" 



^r™(™c)a;fc, 



(11) 



is called symmetry permutation matrix, 
rric is a constant integer called symmetry factor to control 
the maximum length of subsequent positions as center. By 
the way, the component before the subsequent positions and 
consecutive positions to be symmetrized are both created 
randomly. Considering that the symmetry transformation can 
make big change to current solution, it is intrinsically called 
global exploration. For instance, if nic = 0, let choose the 
position 3, then the subsequent position or the center is {0}, 
the consecutive positions {4, 5} with components (3, 1), and 
the function of symmetry transformation is given in Fig. [3] 
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Fig. 3. illustration of the symmetry transformation 

(4) Substitute transformation 

Xk+l = Af''(md)xfc, 



C. Theoretical analysis 

We define the concept of convergence for discrete STA 



\f{xk)-f{x*)\<eyk>N, 



(13) 



(12) 



while, X* is a global optimum of a discrete optimization 
problem, e is a small constant, and A^ is a natural number 
If e > 0, we can say that the algorithm converges to a 
e — suboptimal solution; if e = 0, we can say that the 
algorithm converges to a global minimum. However, if x* is 
a local minimum solution, then we can say that the algorithm 
converges to a e — suboptimal solution and a local minimum 
for e > and e = 0, respectively. 

Theorem 1: the discrete STA can at least converge to a local 
minimum. 

Proof: Let suppose the maximum number of iterations 
(denoted by M) is big enough, considering that the "greedy 
criterion" is used to keep the incumbent best, then there must 
exist a number N < M, when fc > A^, no update of better 
solution will happen. That is to say, f{xk) = f{xN),y k > N, 
where xpf is the solution in the Mh iteration. Denote xn as 
the local minimum solution x*, and then \f{xk) — f{x*)\ = 0. 

■ 

Theorem 2: the discrete STA can converge to a global 
minimum in probability. 

Proof: Let suppose x* = (oi,--- ,a„) is a global 
minimum solution, and xn — (&!,■■ • ,&«,) is the Mh best 
solution. If Xn — x*, according to the "greedy criterion", 
f{xk) = f{xN),y k > N, which means that it converges to 
X* . Otherwise, there must exist a transformation, either swap, 
shift, symmetry or substitute, such that x^+i-^ = (ai, • • • , 6„), 
which means that after li iterations, hi will be changed into 
ai. If /(xjv+ij < /(xat), the x^+h is kept as incumbent 
best for next iteration, else, XN+h is also kept as incumbent 
best for next iteration in probability. Following the similar 

way, after at most I2 + ■ ■ ■ + In iterations, XN+h-{ hi„ will 

be (oi, • • • , a„) in probability. ■ 

D. Parameter selection 

In the state transformations, there are four factors to control 
the intensity between local search and global search. For 
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simplicity and efficiency, the swap, shift and substitute op- 
erators are taken as local search, and the symmetry operator 
is considered as global search; therefore, we consistently make 
ma = 2, TTifc = 1, TTirf = 1 and nic = 0. 

On the other hand, the restore probability pi and the risk 
probability p2 play a significant role in the discrete STA, as 
described by the above theorems. To view the importance of 
the parameters, we arrange a Monte Carlo simulation study. 

Considering the following optimization problem 



where, pc is the penalty coefficient, and p is normally 1 or 2 
(/9 = 1 in this study). Finally, the total cost is 



mill/* 



(14) 



> risk in probability 



t> greedy criterion 



[> restore in probability 



where, /* is created by 

1: r ^0.5,/^r 

repeat 

if f < ri then 

f ^ n 
else if r2 < P2 then 

f '^ n 
end if 
if f* < f then 

f* ^ f 
end if 
if ra < pi then 

end if 
until the maximum number of iterations is met 
here, ri,r2,r3 are uniformly random numbers in (0, 1). 

We test various groups of {pi,p2) for the experiment, in 
which, the maximum number of iterations is le3, and le4 
runs are carried out for each group. The experimental results 
are shown in Table 11] Without loss of generaUty, the group 
(pijP2) — (0.1,0.1) is adopted in this paper for the following 
study due to its good performance and simplicity. 

IV. Implementation of the discrete STA 

The above discrete STA are essentially for unconstrained 
discrete optimization problem. To realize the optimal design 
of water distribution networks, we have to deal with some 
constraints. For the equality constraints on continuity of flow 
and conservation of energy, there exist some hydraulic analysis 
software packages such as EPANET [251, KYPIPE |26|, in 
which the continuity and energy constraints are automatically 
satisfied. Considering that the continuity equations are linear, 
we can first fix some of pipe flows as known to solve the linear 
equations and then substitute them into the energy equations, 
which can reduce the computational complexity of solving 
continuity equations (linear) and energy equations (nonlinear) 
simultaneously. It is not difficult to imagine that the number of 
nonlinear equations equals to that of the simple closed loops or 
required independent paths in a network, and then a Newton- 
Raphson method is used to solve the nonlinear equations. 

For the minimum pressure head constraints, the most com- 
monly used technique is the penalty function method, adding 
a penalty term when the corresponding constraint is violated. 
For example, the following scheme 



J cost Jobj I Jpenal • 



(16) 



NP 
f penal = PC ^ max{0, ifj r 

1=1 



H^Y 



(15) 



A brief description of the steps using discrete STA is given 
in the following 

1) Creat initial Best solution. Generate a group of can- 
didate solutions randomly (the size is the search en- 
forcement, SE) and then select the fittest solution. Let 
Best* = Best and store Best* . 

2) Update the Best. Use swap transformation to generate 
a group of candidate solutions on the basis of Best. If 
the fittest of the candidate solutions is better than Best, 
then accept the fittest solution as Best; otherwise, accept 
the fittest solution as Best in a probability p2- Similar 
procedures are adaptive to shift, symmetry and substitute 
transformations. 

3) Update the Best*. The Best* is updating only when 
Best is better than Best* . 

4) Restore the Best. The Best is restored to Best* in a 
probability pi. 

5) Go back to repeat step 2 until the stopping criterion is 
met. 

Remark 1: We should notice that once a solution is given, 
then the flow in each pipe is determined by solving the nonlin- 
ear equations, and then we can evaluate whether the minimum 
pressure head is satisfied and decide the corresponding penalty 
term to each head pressure constraint. 

V. Case studies 

We investigate the performance of the proposed discrete 
STA by three well-known water distribution networks, namely, 
the Two-Loop network, the Hanoi network and the New 
York network. We first give a detailed study of the Two- 
Loop problem to show that the network system with eight 
unknowns governed by six linear equations and two nonlinear 
equations can be reduced to only two unknowns governed by 
two nonlinear equations. This Two-Loop case is also fully 
trained to study the effect of penalty coefficient and search 
enforcement on the performance of the algorithm. Based on 
the experience gained from the case, the known best solutions 
for the other two networks are also achieved by the algorithm. 



A. two-loop network 

The layout of the two-loop network is given in Fig. |5] There 
are a single reservoir at a 210-m fixed head and eight pipes all 
with 1000-m long. The node data and cost data are given in 
Table III] and Table III and the minimum acceptable pressure 
requirements are all 30-m above the ground level. The Hazen- 
Williams coefficient C is assumed to be 130 for the two-loop 
network. 

In this case study, we give an illustrative procedure of how 
to reduce the complexity of solving the linear and nonlinear 
equations. The flow continuity equations of the two-loop 
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TABLE I 

A Monte Carlo simulation study 



(pi \ P2) 



0.1 



0.3 



0.5 



0.7 



0.9 



0.1 


9.9254e-4± 9.8274e-4' 


0.0010±9.8255e-4 


0.0010 ± 0.0010 


9.9027e-4 ± 9.9498e-4 


9.7423e-4±9.8163e-4 


0.3 


0.0010±9.9359e-4 


0.0010 ± 0.0010 


0.0010 ± 0.0010 


9.8736e-4 ± 9.8159e-4 


9.9254e-4 ± 0.0010 


0.5 


0.0010±9.9755e-4 


9.9249e-4 ±9.7547e-4 


0.0010 ± 0.0010 


0.0010 ± 0.0010 


0.0010 ± 0.0010 


0.7 


0.0010 ± 0.0010 


0.0010 ± 0.0010 


0.0010 ± 0.0010 


9.9456e-4 ± 9.8730e-4 


9.9609e-4 ± 9.9991e-4 


0.9 


0.0010 ± 9.9637e-4 


0.0010 ± 9.9667e-4 


9.8258e-4 ± 9.8261e-4 


0.0010 ± 0.0010 


0.0010 ± 0.0010 



indicates mean it standard deviation 




Fig. 5. two loop network 



TABLE 11 
Node data for the two loop network 



Node 


Demand(m'' 


/h) 


Ground level(m.) 


1 


-1120.0 




210.00 


2 


100.0 




150.00 


3 


100.0 




160.00 


4 


120.0 




155.00 


5 


270.0 




150.00 


6 


330.0 




165.00 


7 


200.0 




160.00 



TABLE 111 
Cost data for the two loop network 



No. 


Diameter (in.) 


Cost ($/m) 


No. 


Diameter (in.) 


Cost($/m) 


1 


1 


2 


8 


12 


50 


2 


2 


5 


9 


14 


60 


3 


3 


8 


10 


16 


90 


4 


4 


11 


11 


18 


130 


5 


6 


16 


12 


20 


170 


6 


8 


23 


13 


22 


300 


7 


10 


32 


14 


24 


550 


Min 


= 2.54 cm 











network are given as follows 

-Q1+Q2 + Q3 + DM2 - 

-Q2 + Q7 + DM3 = 

-Q3 + Qi + Q5 + DA'U = 
-Q7~Qs~Qi + DM5 = 
- Q5 + Qe + DMq = 
-Q6 + Q8 + DM7 ^ 



(17) 



Let Q4, Qq be fixed, then 

Qi == DM2 + DM3 + I?M4 



DM5 + DMf, + DM7 



Q2 = DNh 
Q3 = DNh 
Q5 = DMe 

Qs — Qg ~ 



DM5 + DM7 - Qi 
DMe + Q4 + Q6 

DM7 -Qi-Qe 



Qe 



DM7 



The energy conservation equations can be formulated as 



rsQslQsl 



+ r4Q4\ 
+ r(iQc>\ 



-r7Q7\Q7\ 

-rsQslQs] 



--r2Q2\Q2\ 



(18) 
(19) 



(20) 



and the head loss equations 

H2 = Head-nQilQil'^-' - G2 > i/smin 

H3=H2- r2Q2|Q2r"' - G3 > i?3mi„ 
H4 ~ H2 ~ rsQslQsl'^" — G4 > -ff4min 
H^ ~ Hi — r4Qi\Q4\'^ — G5 > -ffsmin 
Hq= Hi — J'sQslQsT^ — Gq> i?6min 

H7 = Hq~ rgQelQeT" ^ G7 > -ffrmin 

where, Gi{i = 2,- ■ ■ ,7) is the ground level. 

Remark 2: It should be noted that we only need to solve 
the nonlinear system ( [T9| with two unknowns {Qi, Qq). 

For the Two-Loop network, we have to select a diameter 
to each pipe, and for each pipe there are 14 choices. It is 
not difficult to imagine that when choosing a numerical order 
(No.), it corresponds to an exact diameter. That is the reason 
why the discrete STA use the permutation of {1,2, •• • , n} 
as its decision variables and all the intelligent operators are 
operated on a certain permutation. 

Next, we conduct a eremitical study of the the two loop 
network by the proposed discrete STA to investigate the 
influence of the remained parameters, namely, the search 
enforcement {SE) and the penalty coefficient (pc). We set 
SE to be 0.5, 1, 2, 3 and 4 times of the dimension of decision 
variable. Considering that the average cost times the average 
pipe length is 1.0335e5 and the average of the minimum 
pressure heads is 30, the order of magnitude for pc is set at 
Ie4. On this situation, pc is fixed at le4, 2e4, 4e4, 8e4 and le5, 
or increases from le4 to le5 in a linear way. The maximum 
number of iterations is set at 2e2, and a total of 20 runs are 
executed for each group of search enforcement SE and penalty 
coefficient pc. 

As can be seen from Table |IV[ for a fixed SE, the search 
ability is declining as the pc increases, but the feasibility 
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(SE\pc) le4 



2e4 



TABLE IV 
A EMPIRICAL STUDY OF THE TWO LOOP NETWORK 

4e4 8e4 



leS 



Ie4 ■ 



Ie5 



16 



24 



32 



4.2788e5 ± 1.1638e4 

(45%)^ 

4.2015e5 ± 1.2207e4 

(60%) 

4.1826e5 ± 1.4570e4 

(55%) 

4.26 12e5 ± 1.6093e4 

(50%) 

4.2127e5 ± 1.7807e4 
(40%) 



4.3402e5 ± 1.3501e4 

(65%) 

4.2772e5 ± 1.4668e4 

(95%) 

4.2975e5 ± 1.2307e4 

(95%) 

4.2615e5 ± 1.2808e4 

(95%) 

433300 ± 1.6099e4 
(100%) 



4.46 13e5 ± 1.6428e4 

(90%) 

433550 ± 1.3450e4 

(100%) 

434850 ± 1.3880e4 

(100%) 

435400 ± 1.3008e4 

(100%) 

439400 ± 1.0410e4 
(100%) 



445650 ± 1.8930e4 

(100%) 

434900 ± 1.2867e4 

(100%) 

431300 ± 1.17I2e4 

(100%) 

427300 ± 1.5311e4 

(100%) 

434550 ± 1.1432e4 
(100%) 



4.5168e5 ± 1.6885e4 

(95%) 

439000 ± 1.4661e4 

(100%) 

433100 ± I.7669e4 

(100%) 

436100 ± 1.3966e4 

(100%) 

429950 ± 1.2424e4 
(100%) 



4.4757e5±1.6938e4 

(80%) 

4.2901e5 ± 1.3002e4 

(90%) 

430250 ± 1.4400e4 

(100%) 

425600 ± 1.3998e4 

(100%) 

4.37 17e5 ± 1.3846e4 

(85%) 



' indicates the percentage of feasible solutions 



, 


— e— sEsptiri 

— ♦— SBpMel 
— T— SE8p<BM 

— *— 5B8ptl!5 
— *— SBSpcv 


L ■ 



-StlFC2e4 
>SEa|>i:2e4 

-SE2"HjC2*4 
-SE3ipc3e4 




Fig. 6. Iterative curves of best solutions when SE = 8 and pc = 2e4 for the Two-Loop problem, respectively 



rate increases simultaneously with the pc. For a fixed pc, 
the search ability is increasing as the SE increases from 4 
to 8 but declining as the SE increases any more. When the 
pc varies in the iterative process, the performance is not the 
best but much more satisfactory than a constant one to some 
extent. By observation, we can find that setting SE to be 
the dimensionality of the decision variable is a good choice, 
and in this setting environment, pc — 2e4 is a good penalty 
coefficient. Fig. [6] gives the iterative curves of the gained best 
solutions when SE = 8 and pc — 2e4, respectively. It should 
be emphasized that best solutions are all 419,000. 

Remark 3: Under the circumstance, the minimum function 
evaluations to achieve the best known solution is 2048, which 
takes up 0.0001387% of all possible combinations (14*^ = 
1.4758e9). 

Table W\ gives the best solutions gained by various al- 
gorithms, and it can be found that STA can achieve the 
best known solution in this case. It should be noted that 
the same solution was also achieved by GA ["^l, SA pT[ 
and HS |15| with function evaluations at 250,000, 70,000 
and 5,000, respectively. Although the solution in ||6) is even 
better, it should be noted that it brings pipe segments. The 
pressure heads for the Two-Loop network obtained by various 



At node 1, there exists a reservoir with a 100-m fixed head. 



algorithms are given in Table VI 



B. Hanoi network 

The layout of the Hanoi network is given in Fig. U\ There 
are 32 nodes, 34 pipes and 3 loops in this network system. 



The cost data, and pipe and node data are given in Table VII 



and Table |VIII| respectively. The minimum acceptable pressure 
requirements at all nodes are also fixed at 30 m and the Hazen- 
Williams coefficient C is assumed to be 130 as well. 

TABLE VII 
Cost data for the Hanoi network 



No. 


Diameter (in.) 


Cost ($/m) 


1 


12 


45.726 


2 


16 


70.400 


3 


20 


98.387 


4 


24 


129.333 


5 


30 


180.748 


6 


40 


278.280 



From the experience gained from the training of the Two- 
Loop network, the search enforcement SE does not affect the 
performance of the discrete STA explicitly, but the penalty 
coefficient pc plays a significant role in the search ability and 
the solution feasibility, and a good penalty coefficient can be 
evaluated from the order of magnitude the same as the average 
pipe length times the minimum pressure heads. 

For the Hanoi network, the search enforcement SE is set at 
10, and the penalty coefficient pc is fixed at 4e4, or varies from 
4e4 to le5 in a linearly increasing way. The maximum number 
of iterations is set at le3, and a total of 20 runs are executed 
for both fixed and variable pc. Fig. IS] gives the iterative curves 
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TABLE V 

Solutions for the two loop network 



Pipe Alperovits and Shamir 131 Goulter et al. 151 Kessler and Shamir |6| STA (fixed) STA (variable) 



1 

2 

3 

4 

5 
6 

7 
8 

Cost($) 



20 
18 



6 
16 

12 
10 
6 

6 

4 
497,525 



20 
18 
10 

16 
6 

4 
16 
14 
12 
10 
10 
8 
2 
1 
435.015 



18 

12 
10 
16 
3 

2 
16 

14 
12 
10 
10 

8 
3 

2 
417,500 



18 


18 


10 


10 


16 


16 


4 


4 


16 


16 


10 


10 


10 


10 


1 


1 


419,000 


419,000 



TABLE VI 
Pressure Heads for the Two-Loop network 



Node Alperovits and Shamir 1 31 Goulter et al. 5 Kessler and Shamir 161 STA (fixed and variable) 
i_i LJ LJ 



2 


53.96 


54.30 


53.26 


3 


32.32 


33.19 


30.08 


4 


44.97 


44.19 


43.64 


5 


32.31 


32.32 


30.10 


6 


31.19 


31.19 


30.08 


7 


31.57 


31.57 


30.09 



53.27 
30.50 
43.48 
33.85 
30.49 
30.60 




>G>— KD 



G>^^- 



Fig. 7. The Hanoi network 



JOURNAL OF m^X. CLASS FILES, VOL. , NO. , JANUARY 2013 



TABLE VIII 

Pipe and node data for the Hanoi network 



Pipe 


Length (m) 


Pipe 


Length (m) 


Node 


Demand (m? /h) 


Node 


Demand (m? /h) 


1 


100 


18 


800 


1 


-19940 


18 


1345 


2 


1350 


19 


400 


2 


890 


19 


60 


3 


900 


20 


2200 


3 


850 


20 


1275 


4 


1150 


21 


1500 


4 


130 


21 


930 


5 


1450 


22 


500 


5 


725 


22 


485 


6 


450 


23 


2650 


6 


1005 


23 


1045 


7 


850 


24 


1230 


7 


1350 


24 


820 


8 


850 


25 


1300 


8 


550 


25 


170 


9 


800 


26 


850 


9 


525 


26 


900 


10 


950 


27 


300 


10 


525 


27 


370 


11 


1200 


28 


750 


11 


500 


28 


290 


12 


3500 


29 


1500 


12 


560 


29 


360 


13 


800 


30 


2000 


13 


940 


30 


360 


14 


500 


31 


1600 


14 


615 


31 


105 


15 


550 


32 


150 


15 


280 


32 


805 


16 


2730 


33 


860 


16 


310 


- 


- 


17 


1750 


34 


950 


17 


865 


_ 


_ 



of best solutions and changes in 20 runs for the Hanoi problem 
with fixed and variable pc respectively. It is shown that the best 
solution is hit only once by the STA with fixed pc. 

Remark 4: Under the circumstance, the minimum function 
evaluations to achieve the best known solution is 23, 240, 
which takes up 8.1114e-21% of all possible combinations 
(6^4 = 2.8651e26). 



Table IX gives the best solutions gained by various algo- 
rithms, and it can be found that STA with fixed pc can achieve 
the best known solution in this case at the cost of 6.056 million 
dollars, while the solution of STA with variable pc get a 
solution at the cost of 6.065 million dollars. Savic and Walters 
|9_| used the GA to obtain the solution with 1,000,000 function 
evaluations. The solution gained by Zecchin et al. p4) using 
ACO need 100,000 function evaluations. The exactly same 
solution was achieved by SA 1 11 ] and HS 1 15 1 as well, with 
the function evaluations at 53, 000 and 200, 000, respectively. 
The pressure heads for the Hanoi network obtained by various 
algorithms are given in Table IX] 

C New York network 

The layout of the New York network is given in Fig. l9] 
There are 20 nodes, 21 pipes and 1 loop in this network 
system. At node 1, there exists a reservoir with 300-ft fixed 
head. The New York problem is different from other two 
cases, because there already exist pipes in the old system. The 
common objective of this problem is to determine additional 
parallel pipes added to the existing ones to meet increased 
water demands while maintaining the minimum pressure re- 
quirements. The the cost data, pipe and node data are given 



in Table XI and Table XII respectively. The Hazen- Williams 



coefficient C is assumed to be 100 in this case. 

For the New York network, the search enforcement SE is 
also set at 10, and the penalty coefficient pc is fixed at 2e6, 
or varies from le6 to le7 in a linearly increasing way. The 
maximum number of iterations is set at le3, and a total of 




Fig. 9. the New York network 



TABLE XI 
Cost data for the New York network 



No. 


Diameter (in.) 


Cost (S/feet) 


No. 


Diameter (in.) 


Cost (S/feet) 


1 





0.00 


9 


120 


417.0 


2 


36 


93.5 


10 


132 


469.0 


3 


48 


134.0 


11 


144 


522.0 


4 


60 


176.0 


12 


156 


577.0 


5 


72 


221.0 


13 


168 


632.0 


6 


84 


267.0 


14 


180 


689.0 


7 


96 


316.0 


15 


192 


746.0 


8 


108 


365.0 


16 


204 


804.0 
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TABLE IX 

Solutions for the Hanoi network 



Pipe 


Savic and Walters 19 


Zecchin et al. 


14 


Haghighi et al. [l9 


STA (fixed) 


STA (variable) 


1 


40 40 40 


40 


40 


2 


40 40 40 


40 


40 


3 


40 40 40 


40 


40 


4 


40 40 40 


40 


40 


5 


40 40 40 


40 


40 


6 


40 40 40 


40 


40 


7 


40 40 40 


40 


40 


8 


40 40 40 


40 


40 


9 


40 40 30 


40 


30 


10 


30 30 30 


30 


30 


11 


24 24 30 


24 


30 


12 


24 24 24 


24 


24 


13 


20 20 16 


20 


20 


14 


16 12 12 


16 


16 


15 


12 12 12 


12 


12 


16 


12 12 16 


12 


12 


17 


16 20 20 


16 


16 


18 


20 24 24 


20 


24 


19 


20 20 24 


20 


20 


20 


40 40 40 


40 


40 


21 


20 20 20 


20 


20 


22 


12 12 12 


12 


12 


23 


40 40 40 


40 


40 


24 


30 30 30 


30 


30 


25 


30 30 30 


30 


30 


26 


20 20 20 


20 


20 


27 


12 12 12 


12 


12 


28 


12 12 12 


12 


12 


29 


16 16 16 


16 


16 


30 


16 16 12 


12 


12 


31 


12 12 12 


12 


12 


32 


16 12 16 


16 


16 


33 


16 16 20 


16 


16 


34 


20 20 24 


24 


24 


Cost($ millions) 


6.073 6.134 6.190 


6.056 


6.065 



- sift (wlablO 




i3 
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10 
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Fig. 8. Iterative curves of best solutions and changes in 20 runs for the Hanoi problem when pc is fixed and variable, respectively 
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TABLE X 

Pressure Heads for the Hanoi network 
Node Savic and Walters (9) Haghighi et al. | l9| STA (fixed) STA (variable) 



1 


100.00 


100.00 


100.00 


100.00 


2 


97.16 


97.08 


97.17 


97.17 


3 


61.95 


60.82 


61.99 


61.99 


4 


57.21 


56.38 


57.23 


57.34 


5 


51.33 


50.88 


51.31 


51.56 


6 


45.13 


45.13 


45.07 


45.48 


7 


43.68 


43.81 


43.61 


44.06 


8 


41.93 


42.28 


41.85 


42.37 


9 


40.54 


41.09 


40.44 


41.02 


10 


40.34 


37.61 


39.40 


37.01 


11 


38.79 


36.01 


37.85 


35.45 


12 


38.78 


34.83 


34.43 


34.30 


13 


34.58 


30.53 


30.24 


30.10 


14 


36.59 


32.06 


35.49 


33.66 


15 


34.71 


30.96 


33.44 


32.17 


16 


32.08 


31.13 


30.36 


30.53 


17 


33.36 


39.28 


30.51 


33.20 


18 


43.32 


50.04 


44.29 


50.16 


19 


55.54 


57.13 


55.90 


55.37 


20 


50.92 


49.59 


50.89 


50.90 


21 


44.79 


40.04 


41.57 


41.59 


22 


39.63 


34.76 


36.42 


36.44 


23 


44.83 


43.42 


44.73 


44.76 


24 


39.64 


37.73 


39.03 


39.07 


25 


36.38 


34.07 


35.34 


35.40 


26 


32.67 


30.51 


31.44 


31.53 


27 


31.66 


30.32 


30.15 


30.29 


28 


36.48 


38.05 


39.12 


39.15 


29 


32.04 


30.08 


30.21 


30.26 


30 


31.29 


30.58 


30.47 


30.52 


31 


31.81 


30.90 


30.75 


30.80 


32 


32.17 


31.81 


33.20 


33.26 



20 runs are executed for both fixed and variable pc. Fig. 10 
gives the iterative curves of two best solutions and changes 
in 20 runs for the New York problem with fixed and variable 
pc respectively. We can find that the best solution is hit five 
times by the STA with fixed pc and twice by the STA with 
variable pc. 

Remark 5: Under the circumstance, the minimum function 
evaluations to achieve the best known solution is 5200, which 
takes up 2.6883e-19% of all possible combinations (16^^ = 
1.9343e25). 

Table |XIII| gives the best solutions gained by various al- 
gorithms, and it can be found that STA with both fixed and 
variable pc can achieve the best known solution at the cost of 
37.13 million dollars. As a matter of fact, the same solution 
was also gained by GA |9J with the function evaluations at 
1,000,000. The pressure heads for the New York network 



obtained by the discrete STA are given in Table XIV 



VI. Conclusion 

The complexity of the water distribution network comes 
from two aspects, one is the linear and nonlinear equations. 



TABLE XIV 
Pressure Heads for the New York network 

Node STA (fixed and variable) Node STA (fixed and variable) 



1 


300.00 


11 


273.86 


2 


294.33 


12 


275.15 


3 


286.47 


13 


278.12 


4 


284.16 


14 


285.58 


5 


282.13 


15 


293.34 


6 


280.55 


16 


260.16 


7 


278.08 


17 


272.86 


8 


276.51 


18 


261.30 


9 


273.76 


19 


255.21 


10 


273.73 


20 


260.81 



which are commonly handled by a hydraulic solver to ensure 
that the continuity and head loss equations are satisfied auto- 
matically, the other difficulty is that the commercial pipe size 
is discrete, which is proved to be NP-hard. 

In this paper. It is shown that the network system can 
be reduced to the dimensionality of the number of closed 
simple loop or required independent paths, which can reduce 



JOURNAL OF m^X. CLASS FILES, VOL. , NO. , JANUARY 2013 



12 



TABLE XII 

Pipe and node data for the New York network 



Pipe 


Length (feet)** 


Existing Diameters (in.) 


Node 


Demand (feet^/s) ^ 


Minimum Total Head (feet) 


1 


11600 


180 


1 


-2017.5 


300.0 


2 


19800 


180 


2 


92.4 


255.0 


3 


7300 


180 


3 


92.4 


255.0 


4 


8300 


180 


4 


88.2 


255.0 


5 


8600 


180 


5 


88.2 


255.0 


6 


19100 


180 


6 


88.2 


255.0 


7 


9600 


132 


7 


88.2 


255.0 


g 


12500 


132 


8 


88.2 


255.0 


9 


9600 


180 


9 


170.0 


255.0 


10 


11200 


204 


10 


1.0 


255.0 


11 


14500 


204 


11 


170.0 


255.0 


12 


12200 


204 


12 


117.1 


255.0 


13 


24100 


204 


13 


117.1 


255.0 


14 


21100 


204 


14 


92.4 


255.0 


15 


15500 


204 


15 


92.4 


255.0 


16 


26400 


72 


16 


170.0 


260.0 


17 


31200 


72 


17 


57.5 


272.8 


18 


24000 


60 


18 


117.1 


255.0 


19 


14400 


60 


19 


117.1 


255.0 


20 


38400 


60 


20 


170.0 


255.0 


21 


26400 


72 


- 


- 


- 



* 1 feet = 0.3048 m 

' 1 feet^/s = 28.3168 L/s 




- &TA Clixed} 



ItMrtnom 




Fig. 10. Iterative curves of two best solutions and changes in 20 runs for the New York problem when pc is fixed and variable, respectively 



the computational complexity of solving linear and nonlinear 
equations simultaneously to a large extent. 

To overcome the NP-hardness, a new intelligent optimiza- 
tion algorithm named discrete state transition algorithm is 
introduced to find the optimal or suboptimal solution. There 
are four intelligent operators in discrete STA, which are easy to 
understand and to be implemented. The "restore in probability" 
Pi and "risk in probability" p2 strategy in discrete STA is used 
to escape local optimal and increase the probability to capture 
the global optimum. 

At first, a Monte Carlo simulation is studied to investigate 
a good combination of pi and p2, and we find that (^1,^2) — 
(0.1, 0.1) is a good choice. We then focus on a empirical study 
of the Two-Loop network, by training the network, we find that 



the penalty coefficient plays a significant role in the search 
ability and solution feasibility. 

Based on the experience gained from the Two-Loop prob- 
lem, the discrete STA has successfully applied to the Hanoi 
and New York networks, and the results show that the discrete 
STA can achieve the best known solutions with less function 
evaluations. The success of the discrete STA in optimal design 
of water distribution network has demonstrate that the discrete 
STA is a promising alternative in combinational optimization. 
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TABLE XIII 

Solutions for the New York network 



Pipe 


Gessler (2| 


Morgan and Coulter |4| 


Dandy et al. 1 8 


STA (fixed) 


STA (variable) 


1 

















2 

















3 

















4 

















5 

















6 

















7 


100 


144 





108 


108 


8 


100 


144 











9 

















10 

















11 

















12 

















13 

















14 

















15 








120 








16 


100 


96 


84 


96 


96 


17 


100 


96 


96 


96 


96 


18 


80 


84 


84 


96 


84 


19 


60 


60 


72 


72 


72 


20 

















21 


80 


84 


72 


72 


72 


Cost($ millions) 


41.80 


39.20 


38.80 


37.13 


37.13 
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